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We study the exchange and correlation hole of the valence shell of second row atoms using varia- 
tional Monte Carlo techniques, especially correlated estimates, and norm-conserving pseudopoten- 
tials. The well-known scaling of the valence shell provides a tool to probe the behavior of exchange 
and correlation as a functional of the density and thus test models of density functional theory. 
The exchange hole shows an interesting competition between two scaling forms - one caused by 
self-interaction and another that is approximately invariant under particle number, related to the 
known invariance of exchange under uniform scaling to high density and constant particle number. 
The correlation hole shows a scaling trend that is marked by the finite size of the atom relative to 
the radius of the hole. Both trends are well captured in the main by the Perdew-Burke-Ernzerhof 
generalized-gradient approximation model for the exchange-correlation hole and energy. 

PACS numbers: 31.15.ae, 71.15.Mb, 02.70.Ss 



I. INTRODUCTION 

Density functional theory [l], Q , the most widely used 
computational tool for electronic structure calculations, 
is founded upon the knowledge of the existence of a uni- 
versal functional mapping the ground-state density to 
the ground-state energy, with, however, a fundamental 
lack of knowledge how to construct this functional sys- 
tematically. The key problem is describing the effects 
of electron-electron interactions due to Fermi statistics 
and Coulomb repulsion - the exchange and correlation 
energies - in an essentially single-particle description of 
nature. Considerable progress in constructing approxi- 
mate functionals has been made, using a number of varied 
strategies to develop a "Jacobs ladder" hierarchy of mod- 
els of increasing complexity. Two such strategies which 
have proved fruitful are the discovery and implementa- 
tion of scaling laws that describe limiting behavior of the 
universal density functional and the analysis of auxiliary 
expectations, particularly the exchange-correlation hole, 
to provide insight into the role of interelectron correla- 
tions in determining this functional. 

In DFT, scaling laws provide a controlled way of vary- 
ing density, approaching the daunting task of under- 
standing the energy as a functional of the density by first 
tackling the more approachable one of understanding its 
behavior as a function of a scaling variable. Particularly 
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useful is the limit of uniform scaling to high density [H, 0] , 
a process intimately related to the adiabatic connection 
approach to DFT Q , and asymptotically approached by 
the isoelectronic series of atomic ions as nuclear charge Z 
tends to infinity. The properties of exchange and correla- 
tion under this tranformation constrain the possible de- 
pendence of the energy on local-density based variables, 
greatly simplifying the task of functional construction. 
Another fruitful scaling process is that of neutral atoms 
in the Z — > oo limit, in which the charge density tends 
to a well known Thomas-Fermi limit 0, 0] . The gradient 
correction of the latter has been used to diagnose a ma- 
jor limitation in the widely used generalized gradient ap- 
proximation (CCA) of DFT, namely the limited ability 
to tune the GGA to predict accurately both molecular 
and solid-state properties Q, and to motivate effective 
recent remedies for this issue d, [To| • 

Another important tool in the development of DFT 
has been the exchange-correlation (XC) hole. The XC 
hole essentially is the measure of the change in electron 
number density throughout a system given one electron 
observed to be at a given position. The energy of interac- 
tion of an electron with its own exchange-correlation hole 
yields the exchange-correlation energy, the key theoret- 
ical input into DFT. The unexpected degree of success 
of the original local density approximation (LDA) stems 
from the universal properties of the system-averaged XC 
hole (a hole averaged over angle and position in the sys- 
tem) obeyed by the homogeneous electron gas hole @ 
from which the LDA is derived. Input from the behav- 
ior of the XC hole has been important in the develop- 
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ment of effective nonempirical GGA' s [llT [T3| and the 
hybrid DFT-Hartree Fock approach [14J. More recent 
DFT models have not been constructed from XC holes, 
but constructing model XC holes consistent with a given 
functional remains an important analysis tool [l5j |. 

The system-averaged exchange-correlation hole has a 
natural connection to the intracule [l6| or density of 
electron-pairs as a function of their interelectron dis- 
tance. The intracule has been the subject of exten- 
sive study in quantum chemistry, as a source of insight 
into the electron correlation problem. Several classes 
of techniques have been used to study atoms and sim- 
ple molecules, including Hartree Fock (l7l - [T9| . configu- 
ration interaction methods [20Tj25j . and quantum Monte 
Carlo (QMC) [26l - [28j . Among the many applications in- 
clude scaling across isoelectronic series [27| , scaling across 
the periodic table fl7j - Fl9j . decomposition into approxi- 
mate exchange and correlation components [U [2J, [24| 
and shell analysis [l8|, |2(J. Pair densities are of addi- 
tional interest because they can be measured experimen- 
tally [U, [2t| and are the basis of some density functional 
theory approaches [3(1 Hl| . 

The decomposition of the intracule into Fermi, or 
Hartree-Fock hole, and Coulomb hole, incorporating ef- 
fects from correlations, is a close but imperfect equivalent 
to exchange and correlation in DFT. The former uses the 
Hartree-Fock ground-state density as the reference point 
for defining both Fermi and Coulomb correlations rather 
than the exact ground-state density as required in DFT. 
This difference produces a failure of the virial theorem 
at the level of correlation that can be a ten percent ef- 
fect in the case of an atom or a major issue in the case 
of a molecule near dissociation [32j . Secondly, standard 
Coulomb holes account only for the correlation poten- 
tial energy. To obtain the correlation contribution to 
the kinetic energy - the difference between the true and 
noninteracting-system kinetic energies - an adiabatic in- 
tegration of the correlation hole with respect to coupling- 
constant must be performed Calculations of "true" 
exchange-correlation holes, using orbitals that reproduce 
at least the density of one's correlated wavefunction, have 
been done mainly in the QMC approach [33T - [39l | . 

In this paper, we calculate and analyze the exchange 
and correlation holes of the valence shell of second 
row atoms in a pseudopotential model using the vari- 
ational quantum Monte Carlo (VMC) technique. To 
our knowledge, this scaling phenomenon has not been 
studied in the context of density functional theory - al- 
though it was an important ingredient in the earliest at- 
tempts at a pseudopotential description of atomic struc- 
ture Uniform scaling of the radial valence-charge 
density across a row of the periodic table provides a 
convenient way to test the idea of "semilocality" of the 



system-averaged hole by varying the ratio of correlation 
length to atomic radius over a series of otherwise similar 
systems. The second row of the periodic table is easy to 
handle with pseudopotential methods and thus easy to 
isolate strictly valence properties. The valence density 
that results is very close to scale-invariant. The systems 
studied are building blocks of semiconductor materials, 
and probe a density regime, lower than that of the 1st 
row valence shell, typical of many solids for which pseu- 
dopotential simulations are generally used. 

The variational Monte Carlo approach [4l[ makes fea- 
sible the use of explicitly-correlated trial wavefunctions, 
including the electron-electron cusp condition that affects 
the pair-density at zero electron-pair separations. To ob- 
tain separate system-averaged exchange and correlation 
holes, we construct single-particle orbitals that reproduce 
the VMC single-particle density. Exchange holes are then 
calculated numerically from these orbitals. To reduce the 
errors from fluctuations in the pair density due to ran- 
dom sampling, correlated estimates techniques are used. 
These prove important for the measurement of the cor- 
relation hole which is a small fraction of the total pair 
density and thus more affected by noise. The resulting 
exchange and correlation holes are analyzed with respect 
to the scaling of the valence shell density across the row 
and with respect to various density functional models. 

The paper is organized as follows. Section [XT] dis- 
cusses the theoretical underpinnings of the paper - ex- 
change and correlation holes, the GGA approximation, 
and the scaling of the valence shell. Section UTTI describes 
the computation techniques used to generate holes and 
other expectations, Section IIVI presents our results and 
Section [V] our conclusions. All results are expressed in 
hartree atomic units. 



II. THEORY 

A. Expectations of interest 

The exchange-correlation (XC) hole, n xc (r, r + u), is 
defined as the reduction in the ground-state electron den- 
sity from its mean value at some point r + u given the 
observation of an electron at r. It is obtained from a pair 
density fluctuation relationship: 

n(r)n xc (r, r + u) = (r, r + u) — n(r + u)n(r) (1) 

where n(r) is the single-particle density, and 
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is the pair density, measuring the expectation of simulta- 
neously finding electrons at r and r'. Eq. ([l) relates the 
XC hole to the difference between the actual pair den- 
sity and that of the uncorrelated system with the same 
single-particle density. 

The utility of the XC hole in density functional theory 
lies in its relation to the exchange-correlation energy E xc 
through an adiabatic connection 



defined as: 



F _1 

xc ~ 2 



dX 



« 1 

a u —n 

u 



(r,r + u). (3) 



E xc takes into account both the gain in potential en- 
ergy from creating the exchange-correlation hole about 
an electron, and, by means of the integration over the 
coupling constant A, the kinetic energy cost of creating 
the hole as well. The A integral is over a family of systems 
characterized by the same ground-state density but vary- 
ing coupling constant Ae 2 . A A-dependent XC hole n xc is 
defined as an expectation of the corresponding ground- 
state wavefunction. The limits of the integral range from 
a noninteracting system (A = 0), described by the Kohn- 
Sham equations of DFT [42[, to the fully interacting, 
physical system (A = l). 

The A = 1 limit which describes the fully interacting 
Hamiltonian is the focus of this paper. By ignoring the 
integration over coupling constant, we lose the ability 
to calculate the correlation kinetic energy and thus lose 
some of the information contained in the full XC energy. 
However, we keep the ability to assess DFT models for 
this quantity since one can use a scaling relation to con- 
vert the hole integrated over A into that evaluated at 
any specific value of A. This thereby eliminates a tedious 
chore for the quantum Monte Carlo method and allows 
one to explore multiple systems more readily 43| . 

As the Coulomb interaction depends only on the inter- 
particle distance u, and not on the location or angular 
orientation of the hole, it is convenient to define a system- 
and angle-averaged hole, 

(n xc (u)) = Jd 3 r Jdfl u n(r)n xc (r, r + u) (4) 

where d£l u is the solid angle of the pair displacement 
u. This expression may be normalized by the number of 
particles N or number of particle pairs ./V 2 . The system- 
averaged hole contains that part of the XC hole that di- 
rectly affects the determination of E xc , which simplifies 
to 

E xc = [ dXE* c = [ dX [du 2mi{n* c (u)) (5) 
Jo Jo J 

because of the isotropy of the Coulomb interaction. The 
system-averaged hole is closely related to the intracule, 
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The system-averaged hole is the difference between the 
intracule of a quantum system and that of the uncorre- 
lated system with the same single-particle density. 

The XC hole is usefully decomposed in several ways 
to isolate various sources of electron correlation. The ex- 
change (X) hole n x is defined as the exchange-correlation 
hole at zero coupling or n^°. This is the hole associ- 
ated with the Slater determinant wavefunction that re- 
produces the ground-state density of the fully interact- 
ing system and characterizes the inter-electron correla- 
tions due to the Pauli principle. The correlation (C) 
hole is the difference between the exchange and exchange- 
correlation holes and measures the additional many-body 
correlations induced by the Coulomb interaction: 



(n c (u)} = (n xc (u)) - (n x (u)}. 



(7) 



In addition, it is useful to consider the spin- 
decomposition of the XC hole, particularly in spin- 
polarized systems. One may define parallel and anti- 
parallel-spin holes by restricting the sums in Eq. @ to 
particles with the same or opposite spins respectively. 
The parallel spin hole is dominated by the exchange con- 
tribution, since the exchange hole already creates dis- 
tance between electrons of the same-spin, so that further 
effects due to correlation are small. On the other hand, 
the antiparallel spin channel, where the exchange hole is 
zero, contributes the bulk of the correlation hole. 

As it describes the noninteracting system, the spin- 
decomposed exchange-only hole can be written exclu- 
sively in terms of the Kohn-Sham single-particle orbitals 
that define this limit. The system-averaged hole in this 
case is obtained exactly as 



i=l 



V*r(r)C(* + «) 



where 



(8) 



(9) 



and ipia- are Kohn-Sham orbitals for each spin. 

Finally worth noting is that the XC hole obeys the sum 
rule: 



— / 4:iru 2 du(n xc (u)) = -1. 



(10) 



The overall effect of the hole is to remove exactly one 
particle from the measurement of the density about 
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any given electron, essentially removing self-interaction. 
Given their respective definitions, the exchange hole must 
satisfy the same sum rule as exchange-correlation, while 
the correlation hole, being merely a redistribution of the 
N — 1 other electrons around the one in consideration, 
integrates to zero. 



B. Exchange-correlation Hole in Semilocal Density 
Functional Theory 

In a "semilocal" density functional theory, the 
exchange-correlation hole at some point r in space is con- 
structed in terms of various observables defined at that 
point: local spin densities rif (r) and rej (r) in the local 
spin-density (LSD) variant of the LDA [2|,[44(, and adding 
density gradients Vn-^(r)and Vnj (r)in g eneralized gradi- 
ent approximations (GGA's) [l3[ l45rl48j . Kinetic energy 
densities [49U52] ] , and higher-order derivatives of the den- 
sity such as the Laplacian 46, 53, 54] are included in the 
meta-GGA class of theories. 

The LSD exchange-correlation hole at a point r is ob- 
tained from the pair correlation function g xc of the spin- 
polarized homogeneous electron gas (HEG): 

r& SD (r,r + u) = n(r){g* EG [u, r s (r), C(r)] - 1}. (11) 

The pair correlation function is parametrized in terms 
of the Wigner-Seitz radius r s , measuring the average dis- 
tance between electrons, and the spin polarization £, and 
both are evaluated from the local value of the spin den- 
sities: 

\ 1/3 

r«(r) = ( T^TIT (12) 



( 3 
\47rn(r) J 

n(r) 



(13) 



The system-averaged hole and total XC energy may then 
be obtained numerically by applying Eq. Q and Eq. ([5]) 
respectively. 

Among the many variants of the GGA in current 
use, one of particular interest here is that |13j of 
Perdew, Burke and Ernzerhof (PBE), for which models of 
(n x (u)) dHHl] and (n c (u)) [IH have been constructed, 
building upon an accurate HEG hole 56]. Holes within 
the PBE model are designed, under integration, to re- 
produce PBE exchange and correlation energies at any 
value of the density and its gradient within an error of 
about 5%. 

For the exchange hole and energy, the PBE and related 
models introduce into Eq. ([lip a unitless, scale-invariant 
parameter s, defined as 



S (r) 



1 



2k F {v) 



Vu(r) 



n(r) 



(14) 



with kp = (97r/4) 1 / 3 /r s the fermi wavevector. A value of 
s greater than one at a given point indicates the break- 
down of the basic assumption of the LSD that the density 
varies insignificantly on the length scale ~ l/k F of the 
XC hole. For correlation, the PBE employs a second 
inhomogeneity parameter }13j . 



t(r) = [fc F (r)MC)fe s (r)]s(r), 



(15) 



with the local Thomas-Fermi screening vector, k s = 
(4/c_F(r)/7r) 1 / 2 , setting the hole length scale and </>(£), an 
additional scaling factor for spin-polarized systems. 

As the quantity of interest in this paper is the system 
averaged hole rather than the local hole, it is helpful to 
consider system-averages of semilocal DFT parameters. 
A simple definition of the system-averaged Wigner radius 
(r s ) for an atom is [57| 



J dr r 2 n(r) 2 r s (r) 
J dr r 2 n{r) 2 



(16) 



assuming the pair density for zero interparticle separa- 
tion 22j is a reasonable measure of the importance of 



the hole at r to energy expectations, and ignoring the 
anisotropy of the density. System-averaged spin polar- 



and (t 2 



are 



ization (£) and inhomogeneity factors (s 
similarly defined. 

Finally, it is important to note that we are calculati- 
ing the system-averaged correlation hole at full coupling 
constant (A = 1), and with it, the correlation potential 
energy. DFT models for the correlation hole model the 
adiabatically integrated hole and thus the total corre- 
lation energy. It is possible however to construct the 
DFT correlation hole and energy density at a given value 
of coupling constant A from the adiabatically integrated 
version. This can be derived from the scaling properties 
of the hole under uniform scaling of the system [56| . For 
the GGA the following expression for the pair correlation 
function is the result: 



1, * 7 \ - / , , s dg c {r s ,s,(,k F u) 

9 c { r s,s,(,,k F u) = g c {r s ,s,C,k F u) — , 

omr s 

(17) 

with the correlation hole obtained from the expression 



,(r, r + u) = n(r) [g c (r, u) - 1] , 



(18) 



definable for both and n c . The parameters k F u, £, 
and s 2 (but not t 2 ) are invariant under uniform scaling of 
position coordinates Yi (given by — >Ti/X, and n(r) — > 
n(r/A)/A 3 ), and thus held fixed in the derivative. In 
constrast, the exchange hole is invariant under changes 
in scale or coupling constant, and is constructed in terms 
of the scale-independent quantities only. 



4 



REVTeX 4.0 



Cancio, Preprint, 2012 



C. Valence Shell Scaling of Atoms 

A well-known feature of the periodic table is the scal- 
ing of the valence-shell electron density across the 1st or 
2nd row atoms (the 2s and 2p or the 3s and 3p atoms 
respectively). As the number of valence electrons N in- 
creases, the shell radius a(N) shrinks with no change in 
the shape of the distribution: 



n(r;N) 



N 



a(N} 



:n[r/a(N)}- 



(19) 



The scaling behavior strictly occurs for the valence den- 
sity outside the core radius. But with the use of pseu- 
dopotentials to remove the complicated oscillatory be- 
havior of the valence shell inside the core, this scaling 
behavior becomes a global feature of the density - one 
motive for the development of pscudopotentials histori- 
cally [40(. The radius a, usefully defined as the radius of 
the peak in the radial distribution function 4-7rr 2 n(r), is 
a function of valence number N, so that, if the distribu- 
tion were perfectly scaling, the density should reduce to a 
function of TV alone across the row. (We consider neutral 
atoms, with N — Zy, the charge of the ion core.) As a 
result all other expectations of the pseudopotential sys- 
tem ground-state should be reducible to simple functions 
of N, in accordance to the Hohenberg-Kohn theorems [l|. 
Such a picture is limited by the non-scaling behavior of 
the valence density in the ionic core, not a large contribu- 
tion to the probability distribution due to the relatively 
small volume of the core. It is also affected by differences 
between the Mg atom and atoms with 3p orbitals, and 
by the open-shell structure of several of the atoms. 

Important insights into this scaling can be gained by 
the simple heuristic shielding model of the atom devel- 
oped by Slater [4(| • This model assumes a self-consistent 
field felt by an electron in energy shell i given by a 
shielded Coulomb potential (Z—<Ji)/r, with Oi the shield- 
ing charge felt by the electrons of the shell. The orbitals 
for principle quantum number n are nodeless Slater-type 
orbitals r™ -1 exp [{Z — cr^r/n], imposing scaling of the 
density across a row. The effective Bohr radius - the 
peak of the radial probability distribution for the valence 
shell - occurs at 



don 



Z -ay A + BN 1 



(20) 



with ay the valence-shell shielding charge. Slater's 
"back-of-the-envelope" rules for determining shielding 
coefficients yield a shielded effective charge Z — ay that 
is a linear function of N, with coefficients for the valence 
shell of the 2nd row atoms of i? = 0.65 and A=1.55. 

Energetic quantities can thus be constructed as func- 
tions of N from consideration of the scaling form for the 



density [Eq. (fH?)) ] and a(N). The total energy of the 
valence shell in the Slater model is (Z — ay)N/2a(N) 
in hartrees; expressing the shielding charge in terms of 
a(N) yielding a scaling dependence of N/a 2 . The ki- 
netic energy in this naive picture scale in the same way, 
while the external potential due to the ion scales as N 2 /a 
for Zy = N. The Hartree or classic electrostatic energy 
scales as N(N ~ 1) /a provided that self-interaction error 
is removed. With these scaling assumptions, the virial 
theorem relating potential and total energy is satisfiable 
by a (A) taking on the Pade function form of Eq. (I2TJ1) 58] . 

The most interesting energetic quantity for our pur- 
poses is the exchange energy. The heuristic picture for 
scaling of valence-shell energies should be readily ex- 
tended to the exchange energy, as it in fact obeys an 
important universal scaling law Q : 



E x [n a ] = aE x [n]. 



(21) 



where the density n a is defined by the uniform scaling of 
n(r) at constant particle number. 



n a {r) = a 3 n(ar). 



(22) 



To touch base with the current situation, we note that 
particle number N is not fixed as one fills up the valence 
shell, so that the density scales with an additional pref- 
actor N, with an A-dependent length scale a(N) playing 
the role of 1/a. A similar situation is seen in the Thomas 
Fermi scaling of all-electron densities of atoms recently 
revisited in Ref. 8. As with the Hartree energy, and in 
contrast to the all-electron case, self interaction is not 
negligible, varying with particle number and becoming 
critical for the smallest systems. A further complicat- 
ing factor is that as the shell is filled, the overall spin- 
polarization does not stay constant - the process of filling 
the shell is not a uniform scaling of the density itself but 
at best of the radial component of the density. 

To develop scaling model for exchange within DFT, 
we first construct scaling expressions for the system- 
averaged electron-gas parameters (r s ) and (Q defined in 
the previous section. The scaling behavior for (r s ) is that 
of fitting N particles in a box of volume a: 



(r s ) ~ f, = Roa/N 1 / 3 , 



(23) 



with i?o an unknown constant. This equation with 
Eq. (|19p generates similar scaling relations for the in- 
homogeneity parameters s 2 and i 2 , in terms of unknown 
constants Sq and Tq: 



s 2 
t 2 



Sl/N 2/3 
T 2 /aN 1/3 . 



(24) 
(25) 
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A useful estimate of the system-averaged spin- 
polarization (Q is given by 



c 



Nf-N 



N 



(26) 



Note that s 2 does not depend on the scaling parameter 
a, but only on the number of particles N - it is invariant 
under uniform scaling of the density but not, of course, 
under change in N. The polarization £ is also scale in- 
variant. 

A scaling form for the exchange energy can be con- 
structed within the LSD in terms of the scaled parame- 
ters f s and C, incorporating the invariance of exchange of 
the homogeneous electron gas under uniform scaling of 
the density, and a similar spin-density scaling relation- 
ship. The scaled version of the exchange energy within 
the LSD becomes: 



E. r 



where 



<px(0 



N 



0+cf +M) 



~\ 4/3' 



(27) 



(28) 



This is directly obtainable from the definition of the local 
energy-per-particle in the LSD replacing local defi- 
nitions for r s and £ with system-averaged counterparts. 

There is no simple scaling law for correlation to corre- 
spond to that of exchange and the correlation energy is 
difficult to model even for the homogeneous electron gas. 
The analog to the exchange scaling of Eq. ([2T]) docs exist 
for the correlation potential energy [J]: 



Ul\n a ] = aUl /a [n] 



(29) 



tying the process in which the density is scaled uniformly 
by a 3 to that in which the coupling constant is reduced 
at fixed density [6(|. This suggests that the Ar correla- 
tion hole is equivalent to Mg with a electrostatic coupling 
reduced to allow the binding of the full p shell. It does 
not relate the two systems at full coupling. In systems 
where Z is varied at constant particle number, Levy has 
shown [i| that E C (Z) must tend to a constant at large 
Z (that is, uniform scaling to high density). It is possi- 
ble to expect some analogous limiting case for scaling of 
neutral atoms, but this limit is apparently unknown. 



III. SYSTEM AND CALCULATION METHODS: 

A. Hamiltonian 

Our system of interest is described by a many-body 
Hamiltonian for N valence electrons, with a nonlocal 



pseudopotential [6l| to replace the Ne (ls 2 2s 2 2p e ) core: 



N 

E 

i=l 



Xe 2 



N N 

Iz2z2^~- (30) 



The external pseudopotential describing the interaction 
of the valence electrons with the ion core is given by 



V ext (ri) = Vi oc {n) + yWi{n)\lm)(lm 



lr, 



(31) 



including a partially nonlocal term depending upon angu- 
lar momentum projectors \lm). For purposes of compari- 
son to DFT, the Hamiltonian is generalized to consider a 
family of systems characterized by the same ground-state 
density and a variable coupling-constant strength Ae 2 . A 
A-dependent Kohn-Sham potential V^ s is added to the 
external potential to ensure the invariability of the den- 
sity. The range of interest of A varies from zero, describ- 
ing the nonintcracting system, for which the Hamiltonian 
reduces to the Kohn-Sham equation of density functional 
theory, to one, describing the physical system. 



B. Variational Monte Carlo 

For XC hole expectations, we need wavefunctions for 
zero and full coupling respectively. The A = wavefunc- 
tion, i/'Oj is the solution to Eq. (|3"0")) in the abscence of 
electron-electron coupling, and is given by a product of 
Slater determinants of single-particle orbitals. For the 
nonintcracting or Kohn-Sham orbitals we take the output 
of an LSD pseudopotential calculation, with the Kohn- 
Sham potential adjusted to match the spin-densities of 
the A = 1 wavefunction, following the method described 
in Ref . H 

For the A = 1 wavefunction, , we take a variational 
Slater- Jastrow wavefunction of the form 



exp 



(32) 



where D a are Slater determinants composed of orbitals 
from self-consistent LDA orbitals, before adjusting to 
match the VMC density. These are close to, but do 
not exactly match the orbitals of the noninteracting 
wavefunction. Interparticle correlations are described 
through the Jastrow prefactor parametrized by an effec- 
tive pair-potential u. We use a Boys and Handy expan- 
sion of u [62ll63l |. which treats explicitly electron-electron, 
electron-ion, and elcctron-elcctron-ion correlations: 



,(N) 



(r I) r J )= CimnRbirifRbiriTRbira) 71 , 



lmn\l-\-m-\-n<CN 



(33) 
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Here Rb(r) =r/(l + br) and the order of the expansion 
is determined by the factor N = I + m + n. Cusp condi- 
tions [iJ] are used to determine the coefficients C;„ m for 
the linear (N — l) terms; otherwise terms even in TV are 
used and coefficients determined variationally. 

This form of wavefunction is in principle exact for the 
Mg (3s 2 ) spin-singlet system, the analog of the He atom 
for an all-electron model, since there are no spatial nodes 
in the wavefunction and the Jastrow exponential prefac- 
tor contains all possible variables describing the corre- 
lation of three charges. For larger systems, the main 
source of error is likely to be the improper determination 
of the nodal surface of the wavefunction. Multiconfigu- 
rational wavefunctions in which the Jastrow prefactor is 
applied to a linear combination of Slater determinants 
may be useful in this context. Particularly for Al and Si, 
nondynamic correlations in which a (3s 2 ) spin singlet is 
promoted to a (3p 2 ) singlet may be important. At the 
same time, calculations for all-electron systems indicate 
that this kind of effect is important primarily to deter- 
mine the nodal structure between shells, notably, such as 
the Is and 2sp shells of Be (65l - [67| . and should not be as 
important for a single-shell system. 

The variational Monte Carlo (VMC) method [H HI, 
|69| is used to calculate expectations of the trial wavefunc- 
tion and optimize variational parameters. The core of the 
method is to estimate the analytically intractable many- 
body integrals that arise with the use of a Jastrow factor 
by integrating over a randomly selected set of integration 
points R = {ri, r2, • • • , rjy} in the 3iV-dimensional con- 
figuration space. These are sampled with a probability 
proportional to \ipi(R)\ 2 through a a random walk mech- 
anism. Evaluating the nonlocal pseudopotential for an 
integration point requires an additional integration over 
angle for each electron [73, [7lj]> which is done here on 
an 18-point angular grid. Given a set of M such sample 
points, the energy may be estimated by the numerically 
accessible expression 

^ M 

£ ' = M^lHRdBMRi)' (34) 

i 

The variational parameters are determined by the opti- 
mization of the variance of the energy [72j over the sample 
set: 

1 M 

- 2 = mE[^)-^ 2 ( 35 ) 

i 

where E(R) = ^ 1 {R)Hil) T {R). The variance is posi- 
tive definite and approaches zero when the trial wave- 
function globally approaches an eigenfunction, making 
for a robust minimization process with the Levenberg- 
Marquardt algorithm. 



C. Correlated Estimates 

The correlation hole is obtained from differences be- 
tween expectations of the A = wavefunction tpQ and 
that of the fully coupled system tpi, differences which 
may be small enough to make their detection against sta- 
tistical noise difficult, if each expectation were calculated 
independently. Correlated estimate techniques [6^, [73[ , 
specifically developed for calculating arbitrarily small dif- 
ferences in expectations resulting from small changes in 
system parameters or variational parameters, prove to be 
useful here. 

Taking the single-particle density as an example, the 
difference between the radial densities of the two wave- 
functions is obtained from the following expression: 

An(r) = n\(r) — no(r), (36) 

with 

n a (r) = — * 2 , a = 0,1. (37) 

\p \y a {Rk)\ 

k P{Rk) 

The ^-function over the radial distance r can be es- 
timated for a finite sampling set by a histogram 
method (28J. The crucial point for the technique is that 
each term in Eq (1361) is summed over the same set of 
random configurations R = {v\, . . . , rjv}, sampled from 
the probability distribution P(R). By using the same 
random walk for each case, the fluctuations in each eval- 
uation become correlated and can be partially removed 
in taking the difference. In a similar fashion, the system- 
averaged correlation hole can be measured by taking the 
correlated estimate of the difference between the pair 
density expection, or intracule, of the interacting and 
noninteracting wavefunctions. 

This technique does not specify the form of probabil- 
ity distribution P(R) to use in calculating a correlated 
estimate; normally either \ipi\ 2 or |"0o| 2 might be used. 
Either choice can be problematic due to undersampling 
- the instance in which P(R) goes to zero while either tpi 
or ?/>o remains finite. In this situation a rarely sampled re- 
gion of space makes an infinite contribution ~ \ip\ 2 /P to 
an expectation. This can occur if the two wavefunctions 
have different long-range density distributions or nodal 
surfaces, and can ultimately lead to infinite variances in 
the measurement of expectations. 

To eliminate undersampling and optimize the efficiency 
of the calculation, a probability distribution P can be 
chosen consisting of a mixture of the probability func- 
tions for each wavefunction [74] . We use 

P=\(\ij 1 \ 2 -a\^)\+e(\^ 1 \ 2 + a\^ \ 2 )- (38) 
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By mixing both wavefunctions together to form P, it 
becomes impossible for \ipi\ 2 /P ever to become zero - 
should one wavefunction go to zero while the other re- 
main finite, P tends to a nonzero constant. The choice 
of a difference in probabilities emphasizes areas of con- 
figuration space where the difference between the two 
wavefunctions "01 and tpo is largest, and minimizes the 
time spent elsewhere. To obtain a balanced mixture of 
the two states, a is chosen to be roughly the ratio of the 
normalization of the two wavefunctions and may be de- 
termined by measuring the expectation of IV'i/V'ol 2 while 
calculating ground-state energies. The small parameter 
e is chosen to avoid reintroducing undersampling should 
ipi and ipo have the same value. With e^O.l, the method 
reduces the noise in expectations for equal sample sizes 
by a factor of five over non-correlated sampling. 

D. Calculation of exchange hole 

The approach of correlated estimates is well adapted 
for calculation of changes in an expectation between two 
wavefunctions - for our case, the calculation of the corre- 
lation hole. To calculate the expectation for either wave- 
function alone one needs an independent calculation for 
at least one of the wavefunctions. For an atom, it is 
straightforward to calculate the the A = expectation, 
the exchange hole, directly from Eq. ([8]). For an atom 
with only one occupied orbital tp = R(r)Yi m (Q) for a 
given spin, this expression gives a hole proportional to 
the radial component |P(it)| 2 . However for any system 
with two or more orbitals, Eq. ([8]) will involve convolu- 
tions between different orbitals and a more spread-out 
hole. 

This expression is the same up to choice of orbitals as 
the Fermi hole, the system-averaged hole calculated in 
the Hartree-Fock approximation. It is done here numeri- 
cally, using the theoretical formalism of the Hartree-Fock 
intracule calculation of Ref. [HI, Fourier-transforming and 
convolving orbital pairs tpi(r)?p*(r) and summing over all 
possible pairs. It proves important to keep track of the 
angular parts of the wavefunctions so that radial trans- 
forms with high-order spherical Bessel functions are used, 
as well as angular momentum addition techniques. 

IV. RESULTS 

A. Variational calculations and computational 
details 

Table U shows energies and energy variances of vari- 
ational Monte Carlo calculations. We show variational 



Atom 


Wavefunction 


Energy(Error) 


Variance 






(hartree) 


(hartree) 2 


Si 


S 


3.7188(14) 


0.188 


Si 


SJ, N=4 


3.8000(4) 


0.0174 


Si 


SJ, N=8 


3.80422(22) 


0.00943 


Si 


CD, N=8 


3.80524(28) 


0.00942 


Si 


DMC 


3.8065(4) 




Si 


CI 


3.8071 




Mg 


SJ, N=8 


0.84392(4) 


0.000417 


P 


SJ, N=8 


6.52072(25) 


0.0187 


Ar 


SJ, N=8 


21.1922(7) 


0.1039 



TABLE I. Optimized variational energies and variances for 
second row atoms. The symbol SJ is for the Slater-Jastrow 
wavefunction, with the order N indicated; CIJ uses a multi- 
determinant plus Jastrow factor, S is the Slater determinant 
or A = wavefunction, DMC and CI are diffusion Monte Carlo 
and configuration interaction calculations for the same type 
of pseudopotential from Ref. [71Jj . 



energies for the Silicon atom for the Slater determinant 
wavefunction, the Slater-Jastrow correlated wavefunction 
of order N [Eq. (I3"2"j) ]. and a multideterminant- Jastrow 
wavefunction including the low energy 3s 2 to 3p 2 sub- 
stitution into the ground-state determinant [Eq. (|3"3"|) ]. 
These are compared to configuration interaction and dif- 
fusion Monte Carlo results with the same form of pseu- 
dopotential as used here [7l| . Considering the CI results 
as nearly exact, the bulk (92%) of the correlation energy 
has already been achieved for the N — 4 wavefunction. 
The quality of this wavefunction is also indicated by the 
90% reduction in the variance from the noninteracting 
wavefunction. The most accurate wavefunction, includ- 
ing multiple determinants, misses only 2.1% of the corre- 
lation energy, while the most accurate single-determinant 
method, the N — 8 Slater-Jastrow, with 24 variational 
parameters, misses 3.4%. 

The general quality of the optimized wavefunctions 
across the periodic table can be measured by the vari- 
ance of the variational energy, which is zero for the true 
ground-state wavefunction. It is most nearly zero for the 
Mg atom, which is a nodeless wavefunction and in prin- 
ciple exactly treated. The variance for larger systems 
grows faster than the total energy does, but remains rel- 
atively small with the standard deviation in the energy 
for Ar about 1.5% of the total energy. The use of multi- 
configuration wavefunctions proves only to be significant 
for Si and Al. They provide no discernible change in vari- 
ational energy for Mg - consistent with the SJ wavefunc- 
tion being in principle exact for a two-electron system. 
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Atoms with half or more of the p-shell filled lack signif- 
icant low energy configurations and are assumed to be 
described primarily by dynamic correlations. 

Correlation holes are obtained using the correlated es- 
timates technique with ipi taken to be the variation- 
ally optimized N=8 Slater- Jastrow wavefunction and ipo 
the Slater determinant adjusted to give the same single- 
particle density. For 2 x 10 5 samples, the relative sta- 
tistical error in the correlation hole is roughly 1% in the 
physically relevant regime. At long distances, as densities 
and thus Monte-Carlo samples tend to zero, the use of 
correlated estimates eliminate undersampling effectively; 
the worst case statistical errors are about 100% of the 
vanishingly small density differences. For presentation 
in figures, data are processed with a gaussian convolu- 
tion [33| with a width of the order of three histogram 
bins for single-particle data and two bins for the pair 
data. This alters the integrated correlation energy with 
a systematic upwards shift of 1 to 2% but eliminates the 
small residual noise in the curves for better comparison 
between atoms. Energy comparisons use the unsmoothed 
data. 



B. Scaling of the valence density 

Fig. [U shows the scaled radial probability distribution 
Airr 2 an/N for the valence shell of the second row atoms 
from Mg (3s 2 ) to Ar (3s 2 3p 6 ). The scaling parameter 
a is the radius at which the radial probability distribu- 
tion takes on a peak value. The valence density across 
the periodic table clearly reduces to the scaling form of 
Eq. (TT9")) . The scaling is slightly off for atoms Mg, Al, 
and Si with less than half-filled shells but is almost per- 
fect for the other systems. This indicates the insensitivity 
of the radial distribution to the repulsive pseudopotential 
in the ion core region, especially for the larger Z atoms 
for which the core radius shrinks rapidly relative to the 
peak radius a. The relative importance to the density of 
the 3s orbital, with non-zero density inside the core, also 
decreases with increasing Z. 

The atom radius a is plotted in Fig HJa) versus va- 
lence number N; it gradually decreases as the number of 
particles increases, indicating a rapidly increasing den- 
sity. A least squares fit to the Slater shielding model 
[Eq. (|20l) ] with n 2 = 9 for the 2nd row atoms yields an ex- 
cellent fit to a with shielding constants A = 2.350(22) and 
5 = 0.611(4). This corresponds to a 75% effective shield- 
ing of the nucleus from electrons in the 3(s,p) shells by 
each electron in the 2(s,p) shells and 39% shielding by 
each of the other 3(s,p) electrons in the Slater model, 
close to Slater's rule of thumb values of 85% and 35%. 

Also shown in Fig [2] are system averages [Eq. (TT6"]) ] of 
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FIG. 1. (color online) The scaled radial probability distri- 
bution for the valence electrons of the second row atoms Mg 
through Ar, as a function of distance in units of the radius a 
of peak radial density. 



the parameters used in characterizing the semilocal den- 
sity functional theory of the exchange-correlation hole. 
Fig Eta) shows the system-averaged Wigner-Seitz radius 
(r s ). This decreases with N more rapidly than the atom 
radius a, varying from 3 to roughly 0.9 as one goes from 
Mg to Ar. This reflects the two simultaneous effects of 
an increasing number of electrons and a shrinking atomic 
radius as one crosses the row. 

The average polarization (£) shown in Fig [2Kb) ranges 
between zero for the closed shell systems Mg and Ar, 
and a maximum value near 0.6 for P, at half- filling of the 
3p shell, and varies smoothly in between. Although ( is 
defined locally [Eq. (JT3J)], its system average approaches 
the natural global measure, (N t — N I) /N, also shown 
in Fig[2Kb). It might do this better if anisotropic densities 
were used in its system average. 

The inhomogeneity factor s rms — \J (s 2 ) for the ex- 
change hole and t rms = y/ (t 2 ) for the correlation hole are 
shown in Fig (2fc). They show, as expected, somewhat 
different scaling behavior, with s rms decreasing with sys- 
tem size and t rms fairly constant. The values indicate 
systems with a moderate degree of inhomogeneity. As 
one might expect, an isolated atom on average lies out- 
side of the perturbative limit s 2 < 0.3 characteristic of 
solids, but does not quite reach the threshold of severe 
inhomogeneity, s — 1, across the entire system. 

Best-fit global scaled parameters r s , Vs 2 and Vt 2 
[Eqs. (|23|25|26|) ] are also shown in Fig. [2] The fits are 
weighted towards the high-iV end of the row, which shows 
almost perfect scaling behavior. There is excellent agree- 
ment between the actual trend of (r s ) and that predicted 
by scaling, with Rq = 1.592(15). There is also quite good 
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FIG. 2. (color online) Density functional parameters as func- 
tion of valence-electron number N. Part (a) shows the radius 
of peak radial density a and the system-averaged Wigner-Seitz 
radius (r s ), (b) shows the system-averaged spin polarization 
(Q as compared to the fractional difference in occupation 
number (Nf — N±)/N for each spin, (c) shows the system- 
averaged inhomogeneity parameters (t 2 ) and (s 2 ) as compared 
to scaling predictions. 



agreement between the scaling form for (s 2 ) and the ob- 
served value, with a falloff in quality at small N because 
the gradient of the density does not scale as cleanly with 
N as does the density. The value of (t 2 ) shows addi- 
tional behavior that follows the spin polarization. This 
is absent if the system-average is taken over the den- 
sity of antiparallel-spin particle pairs rather than the to- 
tal pair density. Scaling values for these parameters are 
5 = 0.969 and T = 0.856. 




FIG. 3. (color online) Scaled and weighted system-averaged 
exchange hole (a), and correlation hole (b), plotted versus 
scaled distance u/{r 3 ) for Mg, P, and Ar. VMC data are solid, 
the LSD model are dotted and the PBE GGA are dashed 
lines. Curves for P and Mg are shifted downward for clarity; 
in reality they tend to zero for small and large u. 



C. Exchange and correlation holes 

Shown in Fig.[3Ja) are system-averaged exact exchange 
holes evaluated numerically from Eqs. (|SJ) and © and 
DFT models of the same, for three atoms: Mg with a 
filled 3s shell, P with a half-filled 3p shell, and Ar with a 
closed 3p shell. The data for P and Mg have been shifted 
downward slightly for clarity. Each curve is weighted by 
factor of 2iru so that its integral gives the potential energy 
per particle associated with the hole. It thus shows where 
the most significant contributions to the overall energy 
come from. The holes are scaled by a factor (r s ) 3 /N and 
plotted versus the scaled interparticle distance (r s ) to 
reflect the scaling [5(| of a non-polarized exchange hole 
under uniform scaling of the density. The trends for the 
atoms not shown are quite similar. 
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In comparison, the LSD hole [56| is excellent for small 
u, out to around the maximum of the main "dip" that 
contributes the most to the exchange energy. At the same 
time this main feature is too narrow in width. The por- 
tion of the hole missing is spread out into a long range tail 
that reflects the abrupt cutoff in occupation of states at 
the fermi wavevector in a homogeneous electron gas. This 
tail contains a noticeable proportion of the total particle 
sum-rule of the exchange hole (calculated by weighting 
the hole by Ami 2 ) but because of a negligible contribu- 
tion (weighted by 2ttu) to the exchange energy, causes it 
to be underestimated. 

The GGA hole, designed to reproduce the PBE ex- 
change energy [55j , improves upon the LSD by truncating 
the long-range tail and collecting this portion of the LSD 
hole to form a broader hole at its minimum. While failing 
to capture exact details, it mimics quite nicely on average 
what happens in the exact exchange hole, and thus leads 
upon integration to a much improved exchange energy. 
Interestingly, the GGA is slightly less accurate than the 
LSD in the short-u range of the hole, where presumably 
the gradient expansion upon which the GGA is based 
should be most accurate. 

The LSD quite faithfully scales with (r s ) - for example, 
the position of the minimum in the curve is the same for 
all three atoms. The numerical result for the exact hole 
agrees with the LSD for Ar but is shifted to larger u for 
Mg. The GGA only partly follows suit. The reasons for 
this shift will become clear in Section llV E[ 

Figure[3Jb) shows the energy-weighted correlation hole 
for three cases of Mg, P, and Ar. The hole data and 
interparticle distances are scaled and shifted in the same 
fashion as the exchange hole data. The plot shows VMC 
data as a thick solid line, as well as the predictions of the 
LSD and the PBE GGA model [3. 

The correlation holes consist of a short-range region 
where the density of electron pairs is reduced and a re- 
gion at longer distances where it is enhanced; an overall 
sum rule of zero is required. The length scale of the hole 
roughly follows (r s ), increasing as the number of parti- 
cles decreases. The overall hole has a well defined finite 
range, with the density removed at short range collected 
into a noticeable "bump" with a maximum at a distance 
between 1.33 and 2 times that of the valence shell radius 
a. This is intuitively reasonable since there is little phys- 
ical reason to enhance the pair density at interelectron 
distances much larger than the diameter of the atom. 
The shape of the hole varies noticeably from more com- 
pact to more spread out as one moves across the periodic 
table. Likewise, the strength of the correlation hole rel- 
ative to the exchange hole varies considerably, with the 
relative strength of the Mg hole more than twice that of 
Ar. This follows the trend in the homogeneous electron 



gas from highly correlated to uncorrelatcd behavior as r s 
decreases [75[ . Unlike exchange, where the particle sum- 
rule enforces more or less the same size hole in units of 
r s , the zero sum-rule of correlation places little constraint 
on the size of the hole. 

As with exchange, the LSD hole tends to predict the 
short range shape of the hole quite well, with a disagree- 
ment in the on-top (u = 0) hole of 10% hidden by the 
2iru weighting. The hole tends to be too deep and too 
wide. The major disagreement is at long distances: the 
HEG model for correlation includes a long-ranged tail 
that screens out the long-ranged behavior in exchange. 
To compound the effect, the system-average hole at long 
distances is dominated by contributions from the low- 
density asymptotic region of the atom. These in the LSD 
approximation generate unrealistically long tails to the 
system average. The LSD hole for an atom ends up pre- 
dicting that the electron pair density removed at short 
range is redistributed out to infinity at a rate that sur- 
prisingly decays even more slowly than in the HEG. 

The GGA model includes gradient corrections at short 
range, given by a gradient expansion of the HEG model. 
These lift up the correlation hole at short and interme- 
diate distances, creating a markedly better match to the 
VMC results, especially for Ar. The zero sum rule for 
correlation is imposed by a finite-range cutoff of the cor- 
relation hole, which has the added benefit of killing the 
long-range tails of the LSD hole. One therefore finds a 
consistent, systematic improvement on the LSD. How- 
ever the GGA hole dies out too slowly as compared to 
the VMC at long range, so that the positive peak is too 
spread out and contributes less to the energy integral 
than in the VMC case. It is thus easier than in the case 
of exchange to detect the systematic error in the PBE 
model - an overestimate of the size of the correlation 
energy. 



D. Exchange-correlation energy 

Fig. [4] shows correlation, exchange, and exchange- 
correlation energies for our VMC data, and for the LSD 
and PBE density functional models. The VMC data is 
taken from integrating numerically the associated holes 
using Eq. [5] restricted to A = for exchange and A = 1 
for exchange-correlation. The density functional values 
are evaluated directly from the VMC density using the 
energy-per-particle definition conventional for DFT ap- 
plications Hill. The energies are scaled by the exchange 
scaling factor of Eq. (|27p. appropriate, within the LSD 
approximation, for a density which uniformly scales as a 
function of particle number N. This produces a nearly 
constant scaled value for the exchange energy in the LSD, 
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FIG. 4. (color online) Exchange, correlation and exchange- 
correlation potential energy per particle for the valence shell 
of second row atoms, scaled by the factor N <f> x (((,)/ {r s }) as 
described in text. Thick solid line is VMC data; short-dashed 
and long-dashed lines are LSD and GGA predictions. Energy 
in atomic units. Inset: relative differences between the two 
DFT models and the VMC data for the exchange-correlation 
energy. 



indicating the validity of the underlying picture. The 
scaled exchange energy 

E x S ° = -0.390(2) hartree may 
be compared to the homogeneous electron gas value of 
-0.4582 hartree, a reasonable agreement considering the 
arbitrariness inherent in the definitions [Eq. (|16[) ] of (r s ) 
and (C) ■ The spin scaling of the correlation energy is not 
the same as for exchange, as demonstrated by a "bump" 
at half- filling that correlates positively with (£) . 

Fig. 0] demonstrates the dramatic cancellation of er- 
rors in the exchange and correlation components in the 
LSD. The LSD underestimates the effect of exchange, 
having too much of the sum-rule of the hole in its long- 
range tail, and overestimates the resulting screening of 
this long-range tail in correlation. The overall error in 
the exchange-correlation energy, however, is a full order 
of magnitude smaller than that of either exchange or cor- 
relation alone. Having a hole derived from an accurate 
many-body calculation of a true electronic system pays 
off in physically driven error cancellation. The PBE GGA 
implements a consistent correction of these two effects, by 
simultaneously creating a more compact exchange hole 
and cutting off the correlation hole at long range. It rc- 



TABLE II. Mean absolute relative differences (mard) and 
mean absolute difference per electron (mad), in millihartrees, 
between DFT models and numerical data for exchange (X) , 
correlation (C) and exchange-correlation (XC) potential en- 
ergies. 





X 


C 


XC 




LSD 


PBE 


LSD 


PBE 


LSD 


PBE 


mard 
mad 


0.105 
28.4 


0.017 
4.7 


0.657 
26.0 


0.117 
4.6 


0.0076 
2.6 


0.0035 
1.0 



covers the bulk of the errors in the LSD for exchange 
and correlation separately. The averaged difference per 
electron between numerical exchange energies or VMC- 
simulated correlation energies and their DFT counter- 
parts are shown in Table ITT1 The improvement from LSD 
to PBE is an order of magnitude for both exchange and 
correlation but more modest for the two combined. 

The major source of error in our calculation of the ex- 
act exchange energy is from the numerical integration 
of the exchange hole, giving errors typically less than 
one part in 10 5 for the grid used, and is negligible here. 
The VMC results for correlation should in principle be 
a variational upper bound, assuming no numerical er- 
ror in extracting them from pair density data and the 
numerical calculation of exchange. The true correlation 
energy, being lower in energy, should be closer to the 
DFT predictions. To estimate this error in the VMC, 
one can compare the roughly three millihartrcc error of 
our total VMC energy for Si (Table HJ), which we can at- 
tribute to an incomplete treatment of correlation, with 
a 20 millihartree difference between VMC and PBE cor- 
relation energies for Si. On the other hand, the VMC is 
potentially exact for the nodeless Mg valence shell - the 
variational wavefunction used converges to the exact one 
- here the PBE has a 10 millihartree error with respect 
to the VMC, recovering 50% of the error in the LSD. 

The relative difference between VMC and DFT 
exchange-correlation energies is shown in the inset of 
Fig. SJ It is basically on the order of the expected vari- 
ational bias, with a mean signed relative difference of 
less than 0.01%. A expected downwards shift of 5% in 
correlation energy due to variational bias will cause a 
much smaller relative change in the XC energy, given 
that correlation is a small fraction of this energy. The 
shift is 0.5% in the case of Ar; those for smaller N should 
be similar, involving smaller variational biases but larger 
relative correlation energies. Overall, then, the PBE XC 
energy may be roughly 0.5% higher than that of the ex- 
act ground-state, but systematically removing most of 
the LSD error. 
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One may take the exchange scaling analysis one step 
further to analyze the gradient contributions to the ex- 
change and correlation energies. In particular, the GGA 
predicts a multiplicative correction to the LSD exchange 
energy-per-particle of 

^ GA /4 SD ~1 + ^ 2 , (39) 

for small values of s 2 . The PBE data in Fig. 0] shows a 
slight variation from LSD scaling for smaller N or larger 
s 2 which can be fit to Eq. (131)1) . The value of /x thus ob- 
tained is 0.351 for the GGA data as compared to 0.431 for 
the exact calculated exchange, a 20% stronger response 
to inhomogeneity than the GGA prediction. 

E. Scaling trends of exchange and correlation holes 

Under a uniform scaling of the density - and constant 
particle number N - the exchange hole is invariant. As 
our densities scale fairly uniformly, but with N = Z not 
constant, it is informative to see what extent the ex- 
change hole can be reduced to a scaling form. This is 
most easily done by considering the spin-decomposed ex- 
change hole, Eq. ©. To do so, we use an identifiable 
point of each hole - the minimum of the energy- weighted 
hole - to determine a length scale r x for each spin species 
a. The hole is normalized by the number of particles of 
that spin to guarantee a sum rule of -1. Uniform scaling 
then entails the same procedure as for the radial proba- 
bility density: n" x scales to {r x )^n% as distance is scaled 
from u to u/r x . The results are shown in Fig. [SJ 

We do see scaling behavior, but interestingly, two scal- 
ing forms. There is a striking difference between the hole 
formed from a single particle, (in Mg, and the minority 
spin channel of Al through P) and that of two or more 
particles. The two cases that disagree slightly from this 
trend are, quite naturally, the two cases with only two 
electrons in a given spin channel, Al and Si. Otherwise 
the results neatly scale on top of each other. 

This result is not that surprising if we consider the 
form of the exchange hole in Eq. For one parti- 

cle in a given spin, the form reduces to a convolution of 
the spin density n a (r) = \tp^ s (r)\ 2 and leads to a rela- 
tively compact function. The hole merely removes the 
self-interaction error of the Hartree approximation, and 
in a sense is not a true exchange hole. The N a > 1 case 
also includes convolutions of the overlap of two different 
orbitals tp* a (r)tpj^(r), that describe the exchange of two 
electrons. Such overlap terms are naturally more spread 
out than a single-orbital probability and create a more 
slowly decaying hole. Note that the exchange scaling law 
[Eq. (1211) ] requires that both forms scale uniformly with 
an isolectronic (fixed N) uniform scaling of the density, 




0.0 1.0 2.0 3.0 4.0 5.0 6.0 



u/l ax 

FIG. 5. (color online) Scaled minority-spin (thinner lines) 
and majority-spin (thicker lines) exchange holes for second 
row pseudopotential atoms plotted versus a spin-dependent, 
scaled interparticle distance u/r x . Each hole is scaled by 
(r x ) 3 and weighted by 2n times the scaled interparticle dis- 
tance. The scaling length is determined by the distance 
at which the weighted hole has its minimum value. 



but with significantly different asymptotic forms because 
of their different origins. Finally it is interesting to note 
how quickly the transition from the self-interaction dom- 
inated to an exchange dominated hole occurs - the large 
number limit is essentially reached starting at N a = 2. 

Fig. [5] shows scaling lengths for the majority (up) spin 
and minority (down) spin exchange holes, r x and r x , 
in units of the valence-shell radius a for each atom in 
the second row. In addition, the average Wigner-Seitz 
radius (r s ) is plotted and the equivalent spin-dependent 
radii (rj) = (?' s )(l + erC) -1 / 3 , proportional to the natu- 
ral length-scale (kp) -1 of the spin-decomposed LSD ex- 
change hole. These are scaled by a factor of 0.755 so that 
they match r x for Ar. The comparison between the ac- 
tual length-scales rx and the LSD equivalent shows the 
separation of the single-particle and multi-particle cases 
seen in Fig. [5] For spin occupation N a > 1, the rx values 
are well predicted by LSD theory. For N a = 1 the hole 
scales as a, and is notably larger than the LSD predic- 
tion. In this case, the length scale is set simply by the 
width a of the single orbital occupied. 

Unlike exchange, the universal scaling law for correla- 
tion [Eq. (I29[) ] is nontrivial, and the correlation energy 
and hole are nontrivial to model even for the homoge- 
neous electron gas. Nevertheless, Fig. |3Jb) shows that 
correlation holes for the second row atoms are qualita- 
tively similar and it is instructive to scale the holes to 
highlight the trends that occur as one crosses the row. 

To do this we use an empirical scaling relation 
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FIG. 6. (color online) Various scaling lengths discussed in 
the paper. Shown versus valence electron number are: the po- 
sition of maximum depth (rx) in the spin- up (red) and spin- 
down (blue) energy-weighted exchange holes, the same (re) 
for antiparallel (green) and parallel spin correlation holes, all 
scaled by the effective valence bohr radius a. System-averaged 
Wigner-Seitz radius (r s ) for the total density (black) and each 
spin density (dotted lines) are also shown, scaled by a factor 
of c = 0.755 to aid in comparison. 
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that matches the holes as closely as we find possible. 
It is again helpful to spin-decompose the holes, into 
antiparallel-spin (A) and parallel-spin (P) channels. A 
scaled form n c (x) of the hole in either channel may be 
constructed as 

2TTu(n^ p {u)) = 2Tr(u/r£ P )N A , P nf> p (u/r£ P ), (40) 

where Na and Np are the number of electron pairs of 
either spin channel and and are scaling lengths. 
These are again chosen as the distance at which the 
energy- weighted holes take their minimum value. This 
optimizes the match between holes at short interparti- 
cle distances at the expense of that at longer ones. The 
results are shown as a function of scaled interparticle dis- 
tance in Fig. [7] scaling lengths are shown in Fig. [6] 

For the case of antiparallel-spin electron pairs 
[Fig-E^a)] , the results at short distances match up closely. 
The minimum of the holes shows some shell structure, 
with the closed spin-shell atoms Mg, P and Ar with the 
deepest holes. At long range, there is a systematic trend 
as the 3p-shell is filled, going from a compact hole with 
a sharp positive peak, to a relatively wider hole with a 
positive peak spread out over a large range of distance 
relative to r^. The trend seems to be one of gradual re- 
duction of finite-size effects at long-range, with the shape 
of the hole trending to that of the homogeneous electron 
gas, with its infinite-ranged hole. 



x=u/r p 

FIG. 7. (color online) Scaled correlation holes for second 
row pseudopotential atoms. Part (a) shows the correlation 
hole at full coupling for antiparallel-spin electrons, weighted 
by 2ir times the interelectron distance u. The holes are scaled 
by the number of antiparallel electron pairs Na and plotted 
versus the scaled electron distance x = u/r A , where r A is the 
distance at which each weighted hole has a minimum. Part 
(b) shows the analog for the parallel-spin correlation hole, in 
terms of parallel-spin pair number, Np and scaling radius r^. 
Each hole is plotted in atomic units. 



The parallel-spin holes are shown in Fig. EJb). That 
for Mg is trivially zero since there are no same-spin pairs. 
As shown in Fig. [HI r p is almost exactly twice that of 
across the second row so that each subplot of Fig.[7]shows 
the same physical range in distance for each atom despite 
the different abscissas. The hole-per-pair is significantly 
smaller than in the antiparallel spin channel and is con- 
centrated at longer interparticle distance. Both effects 
are caused by the exchange hole which removes most 
of the probability of finding particles of the same spin 
at short range. As one goes across the shell, the ten- 
dency is for gradually deeper and longer-ranged holes. 
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Finally, given the equivalence r^ — 2r^, one finds that 
the scaled antiparallel and parallel holes-per-pair for each 
atom match fairly closely at large distances (tq > 5) . At 
distances longer than the effective exchange hole radius, 
it seems that electrons no longer detect each other's spins. 

All-electron spherically-averaged Coulomb hole data 
for nearly spherical molecules and atoms show quali- 
tatively very similar results to ours, after taking into 
account the Airr 2 weig hting typically used in the liter- 
ature [13, E3, [If. The hole at small interparti- 
cle distances typically shows the evidence of shell and 
molecular structure, but significantly reduced by the ef- 
fects of system averaging, particularly in the large-Z 
limit. The qualitative shape of the peak at large inter- 
particle distances seems to be a universal feature of the 
Coulomb/correlation hole and closely matches the be- 
havior seen in Fi gs. Ef b) and [7] Shell analysis of the hole 
for small atoms [2(| tells us that this feature is caused 
by electrons in the valence shell and is thus apropos for 
comparison to our data. Using the point of crossover 
r cr from negative hole to positive peak as a point of ref- 
erence, we find that our scaling model provides a fairly 
good prediction for other data. A fit of this point to our 
data yields r cr = 1.5*Roa/N 1 ^ 3 for Si; using this formula 
and the naive Slater model for a we predict Coulomb hole 
crossover points for the C isolectronic series that range 
from 11.2 a.u. for C to 8.3 a.u. for Ne +4 , which are 
indistinguishable from those of the Coulomb hole plots 
of VMC data for the series [27| • This suggests that our 
data may be useful to analyze the long-range trends of 
correlation holes of atoms and small molecules. 

A few extra notes are in order - first of all, the point of 
comparison between holes found to be most successful is 
the correlation hole per pair and not per particle, as with 
exchange. A particle number for the antiparallel-spin 
hole cannot be unambiguously defined, and we find no 
scheme for a per-particle hole that provides as consistent 
a scaling fit for the series. At the same time, the scaled 
hole used here is not scale-invariant, in which case the 
hole, having units of volume, should scale as l/r^ cale . It 
is therefore not unitless, but rather varies as l/ctp given 
the choice of atomic units. The situation is reminiscent 
of that of the slowly varying electron gas, in which the 
inhomogeneity in the exchange hole can be described by 
the scale-invariant parameter s, but the correlation hole 
depends on the non-scale-invariant t ~ s/y/r^. However, 
it is not impossible to scale the correlation hole with a 
uniform scaling hypothesis of the form: 



(n c {u)) 



F(N t ,NJ 
(r 



scale ) 



n c (u/r s 



(41) 



scaling behavior is assumed to lie in the arbitrary ampli- 
tude F. The two approaches may be considered equiva- 
lent since the normalization used in Eq. (|40l) is implicitly 
a function of Nf and N±. 

It is also worth noting that length-scales for correlation 
and exchange diminish as a proportion of atom radius a 
in going from Mg to Ar (Fig. [HI . Finite size effects in ex- 
change and correlation become important as (r s )/a ap- 
proaches unity - that is, the length-scale of the XC hole 
approaches the system size of the atom. We thus expect, 
and find, the largest errors for local DFT's for Mg, for 
which (r s )/a is largest, and the smallest for Ar. Nev- 
ertheless, the GGA parameter t 2 used to estimate this 
inhomogeneity error has a system average that increases 
as one proceeds down the row. In the slowly varying 
electron gas, the higher the electron density, the more 
sensitive correlation is to inhomogeneity, but in atoms, 
the higher the density, the less effect the finite size of the 
atom has on correlation. This misidentification is a natu- 
ral limitation of using a semilocal parameter to measure 
the effect of inhomogeneity, which for atoms as essen- 
tially zero-dimensional objects must necessarily depend 
on global features. Interestingly, however, the PBE cor- 
relation hole does take into account one global measure, 
the zero particle-sum-rule of the hole. For very inhomo- 
geneous systems, this constrains the hole to react more 
strongly at larger r s to a given value of t 2 . Thus it fol- 
lows to a fair degree the observed finite-size effects, more 
so than one might have expected. 



F. Pair correlation function 

One can gain more insight into the behavior of corre- 
lation with atomic number, and in particular the role of 
finite size effects, by calculating a pair correlation func- 
tion g c , defined as the ratio of the pair density of the 
fully correlated (A = 1) system to that of the equivalent 
noninteracting system (A = 0): 



,(2) 



9c(u) 



i(«)> 



o(«)> 



(42) 



for an arbitrary function F(N^, iVjJ of the number of va- 
lence electrons of each spin. In this case, the absence of 



This maps the fractional change, because of Coulomb 
correlation, in the expected number of particle pairs as a 
function of distance. As the noninteracting-system pair 
density already incorporates exchange, this measures a 
purely correlation effect, and is somewhat different from 
the normal definition of g measured relative to the pair 
density of independent particles. The value at zero sep- 
aration, g c (0), measures the on-top hole for antiparallel 
spin correlation. Because of fermi statistics, the parallel- 
spin channel contributes zero to both numerator and de- 
nominator of the on-top value of Eq. (|42|) . 
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u (a.u.) 



FIG. 8. The correlation contribution to the pair correlation 
function, g c , for the valence shell of the second row atoms 
plotted versus interparticle distance it. The curve for Ar and 
Mg are labelled at u = 0. Also shown (circles) is the system- 
averaged on-top or u = pair-correlation function for the LSD. 



We show the VMC values for g c (u) for the second 
row atoms in Fig. [8] As a comparison, the equiva- 
lent system-averaged on-top pair correlation function for 
the LSD, (gc SD (0)), calculated by the same system- 
averaging technique as Eq. (|16[) . is shown as circles at 
u = 0. At short range the pattern of the VMC data is 
similar to the correlation hole in the HEG 75] except 
for deviations, due to small-number statistics, from an 
expected cusp [76| at very short interparticle distances. 
Given the noise in the VMC in this region, the LSD and 
VMC values for the on-top hole are in reasonable agree- 
ment. At long range one sees a marked enhancement of 
the distribution of particle pairs relative to the uncor- 
related pair density. The long-range asymptotic value 
of the pair distribution function in the HEG is one, in- 
dicating a vanishing difference between correlated and 
uncorrelated distributions. In atoms, this is not a re- 
quirement, and in fact the enhancement of pairs at long 
range is surprisingly large for Mg. The distances in this 
case are bigger than the distance from the peak of the 
valence density on one side of the atom to that on the 
other, about 1.7 a or roughly 3.5 a.u. for Mg. There is a 
large fractional enhancement of the relatively small den- 
sity of pairs separated by several atomic radii, in contrast 
to the modest enhancement of the long-range pair den- 
sity for more localized holes. This is consistent with an 
"in-out" correlation where if one electron is found on the 
inside of the valence shell the other is favored to be found 
on the outside edge of the shell. 



V. DISCUSSION 

The valence shell of first or second row atoms, within 
the pseudopotential approach, is an example of uniform 
scaling that has been known from the earliest stages of 
atomic physics. The form of scaling is related to that of 
the Thomas-Fermi scaling of the all-electron density for 
large atoms, in that the net charge of the system is kept 
fixed as the system is scaled, and so scaling parameters 
such as (r s ) depend upon the number of electrons N. 
It differs in that the scaling parameter a does not have 
a simple power-law behavior with N, but essentially a 
Pade-like dependence, as noted by Slater, because of the 
role of self-interaction. At the same time, the scaling of 
a single shell is a form of uniform scaling to high density, 
and might be a useful complement to the standard exam- 
ple of isoelectronic scaling. What happens in the current 
case, in the limit of a very large-degeneracy shell, for 
which N ^> 1 is achievable? 

Intriguingly, the exchange hole has not one, but two, 
scaling forms - with the single-orbital hole fundamen- 
tally different from those that are constructed with two 
or more orbitals. Self-interaction effects spoil the invari- 
ance of the exchange hole under uniform scaling with con- 
stant net charge that would hold with constant particle 
number. But at the same time, self-interaction effects 
die out astonishingly rapidly, with systems of three or 
more electrons already showing large-iV scaling behav- 
ior. Correlation holes do not scale uniformly and trends 
across the second row cannot be unambiguously mod- 
eled. Nevertheless, a clear trend with the scaling of the 
density does occur: the transition from more correlated, 
low density systems in which finite size effects dominate 
the correlation peak at large interparticle distances, to 
higher density ones in which the system-averaged hole ap- 
proaches that of the HEG, and the GGA approximation 
is very good. The natural parameter to characterize this 
trend, the ratio (r s )/a of average hole radius to atomic 
radius, can not be modeled in a semilocal DFT, but the 
use of a different constraint - a cutoff based on the corre- 
lation hole sum rule - does a reasonable job of following 
it. The qualitative behavior of the positive correlation 
peak is seen in all-electron calculations and a crude esti- 
mate indicates that scaling results for this feature should 
apply to atoms and spherical molecules in general. 

It is of interest to analyze the sources of error in the 
PBE correlation hole by decomposing the system aver- 
aging. The PBE hole near the valence peak determines 
the overall shape of the system-averaged hole but, since 
gradients are small in this region, has an unphysical long- 
range tail equivalent to that of the LSD hole. Holes in the 
pseudopotential core and far from the atom deviate dra- 
matically from the system average but cut off just at the 
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right length scale. They end up contributing the most to 
the positive peak of the system-average that, physically, 
is caused by finite size of the atom. It is not the local XC 
hole but its system-average that the PBE is capturing, 
as it was designed to do. 

Our data thus show that for the type of system studied 
here, the semilocal PBE GGA model works essentially as 
advertised. The significant defects of the XC hole in the 
LSD approach are more or less fixed, especially for the 
crucial aspect of the hole - the integral that produces the 
exchange and correlation energy. The reason for this suc- 
cess is likely related to the simplicity of the single-shell 
structure studied, with the main physics being scaling 
behavior similar to that which the PBE was built to rep- 
resent. The main source of error, the poor treatment of 
finite-size effects, occurs mainly in the long-range tail of 
the GGA XC hole which does not contribute much to the 
total energy. 

The PBE does somewhat underestimate the gradient 
correction parameter \x [Eq. Q39p] needed for valence-shell 
exchange energies. Recent work on modifying the PBE 
for solids [43, |48| indicates that use of a value for /i 
half that of the PBE, but consistent with the gradient 
expansion for a slowly- varying electron gas, leads to im- 
proved lattice constants and bulk moduli (if poorer cohe- 
sive energies). The PBE choice, is instead best suited for 
predicting total energies of atoms and binding energies 
of molecules. Our work thus emphasizes the incompat- 
ibility between GGA's designed for molecules and those 
for solids. The large value of \i needed for total atomic 
energies has been attributed [8| to using the gradient cor- 
rection to account for exchange energy corrections caused 



by the cusp in electron density near the nucleus. This is 
should not be true in the present case since the nucleus 
has been replaced by a smooth pseudopotential. It seems 
that another mechanism is at play here, quite possibly 
self interaction. The PBE models systems with large self 
interaction error, like the two-electron valence shell of 
Mg, very well. But this is perhaps at a cost of overcor- 
rccting in situations like the slowly varying electron gas 
where inhomogeneity and not self-interaction is the pre- 
dominant issue. Another consideration is the boundary 
condition difference between finite and infinite systems, 
with the former needing a stronger correction to produce 
a more sharply defined long-range cutoff to the hole. 

Our work also points to the difficulty of imposing a 
self-interaction correction to the GGA. To the extent 
that the GGA may be correcting LSD error caused by 
self-interaction and not by the imperfect treatment of in- 
homogeneity, applying a self-interaction correction to the 
GGA would lead to correcting the same problem twice. 
Any self-interaction correction based on a GGA would 
require a remarkable degree of cancellation of errors be- 
tween the correction for exchange and that for correlation 
to improve total energies for the systems studied here. 
The exploration of self-interaction corrected GGA mod- 
els that could have this level of error cancellation might 
thus be of interest. 
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